import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.arima_model import ARIMAResults
import warnings
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.iolib.smpickle import load_pickle
import statsmodels.api as sm

plt.rcParams['font.sans-serif'] = ['SimHei']  # 定义使其正常显示中文字体黑体
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示表示负号

warnings.filterwarnings("ignore")

data = pd.read_csv('AirPassengers.csv', index_col=u'Month', header=0)
print(data.head())

# 画出时序图
data.plot()
plt.show()

# 画出自相关性图
plot_acf(data)
plt.show()

# 平稳性检测
# adf : float
#     The test statistic.
# pvalue : float
#     MacKinnon's approximate p-value based on MacKinnon (1994, 2010).
# usedlag : int
#     The number of lags used.
# nobs : int
#     The number of observations used for the ADF regression and calculation
#     of the critical values.
# critical values : dict
#     Critical values for the test statistic at the 1 %, 5 %, and 10 %
#     levels. Based on MacKinnon (2010).
# icbest : float
#     The maximized information criterion if autolag is not None.
# 单位检测统计量对应的p 值显著大于 0.05 ，说明序列可以判定为 非平稳序列
print('原始序列的检验结果为：', adfuller(data[u'#Passengers']))

# # 对数据进行差分后得到 自相关图和 偏相关图
# D_data = data.diff().dropna()
# D_data.columns = [u'#Passengers差分']
# D_data.plot()  # 画出差分后的时序图
# plt.show()
# plot_acf(D_data)  # 画出自相关图
# plt.show()
# plot_pacf(D_data)  # 画出偏相关图
# plt.show()
# # 一阶差分后的序列的时序图在均值附近比较平稳的波动， 自相关性有很强的短期相关性，
# # 单位根检验 p值小于 0.05 ，所以说一阶差分后的序列是平稳序列
# print(u'差分序列的ADF 检验结果为： ', adfuller(D_data[u'#Passengers差分']))  # 平稳性检验
# # 对一阶差分后的序列做白噪声检验
# print(u'差分序列的白噪声检验结果：')  # 返回统计量和 p 值 , lags为检验的延迟数
# print(acorr_ljungbox(D_data, lags=20))
# # 差分序列的白噪声检验结果： (array([*]), array([*])) p值为第二项， 远小于 0.05
# # 如果pvalue值大于0.05，就说明我们无法拒绝原假设（该序列是白噪声序列），这个时间序列的确就是白噪声序列,无分析价值


data["Passengers_diff1"] = data[u'#Passengers'].diff(1).dropna()
data["Passengers_diff2"] = data["Passengers_diff1"].diff(1).dropna()
data1 = data.loc[:, ["#Passengers", "Passengers_diff1", "Passengers_diff2"]]
data1.plot(subplots=True, figsize=(18, 12), title="差分图")
plot_acf(data["Passengers_diff1"].dropna(), title="Autocorrelation Passengers_diff1")  # 画出自相关图
plot_pacf(data["Passengers_diff1"].dropna(), title="Partial Autocorrelation Passengers_diff1")  # 画出偏相关图
plot_acf(data["Passengers_diff2"].dropna(), title="Autocorrelation Passengers_diff2")  # 画出自相关图
plot_pacf(data["Passengers_diff2"].dropna(), title="Partial Autocorrelation Passengers_diff2")  # 画出偏相关图
plt.show()
# 单位根检验 p值小于 0.05 ，所以说一阶差分后的序列是平稳序列
print(u'1阶差分序列的ADF 检验结果为： ', adfuller(data[u'Passengers_diff1'].dropna()))  # 平稳性检验

# 对一阶差分后的序列做白噪声检验
# 差分序列的白噪声检验结果： (array([*]), array([*])) p值为第二项， 远小于 0.05
# 如果pvalue值大于0.05，就说明我们无法拒绝原假设（该序列是白噪声序列），这个时间序列的确就是白噪声序列,无分析价值
print(u'1阶差分序列的白噪声检验结果：')  # 返回统计量和 p 值 , lags为检验的延迟数
print(acorr_ljungbox(data[u'Passengers_diff1'].dropna(), lags=10))

print(u'2阶差分序列的ADF 检验结果为： ', adfuller(data[u'Passengers_diff2'].dropna()))
print(u'2阶差分序列的白噪声检验结果：')
print(acorr_ljungbox(data[u'Passengers_diff2'].dropna(), lags=10))

# 对模型进行定阶
D_data = data[u'Passengers_diff2'].dropna()
pmax = int(len(D_data) / 10)  # 一般阶数不超过 length /10
qmax = int(len(D_data) / 10)
bic_matrix = []
for p in range(pmax + 1):
    temp = []
    for q in range(qmax + 1):
        try:
            mod = sm.tsa.arima.ARIMA(D_data, order=(p, 2, q))
            res = mod.fit()
            # res.summary()
            value = res.bic
            temp.append(value)
        except Exception as e:
            print(e)
            temp.append(None)
    bic_matrix.append(temp)

bic_matrix = pd.DataFrame(bic_matrix)  # 将其转换成Dataframe 数据结构
print(bic_matrix)

# bic_matrix.stack(), index为MultiIndex，【行，列】
p, q = bic_matrix.stack().idxmin()  # 先使用stack 展平， 然后使用 idxmin 找出最小值的位置
# p, q = 2, 4
print(u'BIC 最小的p值 和 q 值：%s,%s' % (p, q))  # BIC 最小的p值 和 q 值：2,4

# 所以可以建立ARIMA 模型
mod = sm.tsa.arima.ARIMA(D_data, order=(p, 2, q))
res = mod.fit()
# 保存模型
res.save('model.pkl')
print(res.forecast(5))

# 模型加载
loaded = load_pickle('model.pkl')
# 预测未来五个单位
predictions = loaded.forecast(5)
print(predictions)
